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We analyze the spectrum of axions radiated from collapse of domain walls, which have received 
less attention in the literature. The evolution of topological defects related to the axion models 
is investigated by performing field-theoretic lattice simulations. We simulate the whole process of 
evolution of the defects, including the formation of global strings, the formation of domain walls 
and the annihilation of the defects due to the tension of walls. The spectrum of radiated axions has 
a peak at the low frequency, which implies that axions produced by the collapse of domain walls are 
not highly relativistic. We revisit the relic abundance of cold dark matter axions and find that the 
contribution from the decay of defects can be comparable with the contribution from strings. This 
result leads to a severer upper bound on the axion decay constant. 
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I. INTRODUCTION 



The standard model of particle physics has been well established, except several unresolved problems. Among 
them, the strong CP problem in quantum chromodynamics (QCD) remains as a mystery for a long time. The most 
attractive solution is the celebrated Peccei-Quinn (PQ) mechanism [l|, which introduces an global U(l) symmetry 
that has to be spontaneously broken at some high energy scale. This spontaneous breaking of the global symmetry 
predicts an existence of a (pseudo) Nambu-Goldstone boson, called the axion 0. Soon after the proposal, it was 
found that the prototype model of the axion conflicts with experiments 0], and it was ruled out. However, it was 
argued that models with higher symmetry breaking scale denoted as F a (the axion decay constant) can still avoid the 
experimental constraints (HQ ■ The essential point is that the couplings between axions and other fields are suppressed 
by a large factor of the symmetry breaking scale ~ l/F a . These models are called "invisible axions" because of their 
smallness of coupling with matter. 

This invisibleness leads to a cosmological consequence. It turns out that almost stable coherently oscillating axion 
fields play a role of the cold dark matter filled in the universe The present density of the cold dark matter is 
measured to be Ocdm ~ 0.23 0- Requiring that the present abundance of axions must not exceed the observed 
value, one obtains the upper bound for the symmetry breaking scale F a < 10 12 GeV. In other words, axions are good 
candidate of the dark matter if F a is as large as 10 12 GeV. 

Furthermore, the astrophysical observations give a lower bound on F a . For example, from the helium- burning 
lifetimes of horizontal branch stars in the globular clusters one can obtain a bound F a > 0(1) x 10 7 GeV |[. Also, the 
cooling time of white-dwarfs @ gives a bound F a > 0(1) x 10 8 GcV for Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) 
model |5|. The most stringent bound comes from the energy loss rate of the supernova 1987A [3, which gives 
F a > 10 9 GeV both for Kim-Shifman-Vainshtein-Zakharov (KSVZ) model and for DFSZ model. Combined with 
the cosmological bound described above, we can constrain the axion models into the parameter region with F a ~ 10 9 - 
10 12 GeV. This is called "the classical axion window". 

It was argued that this axion window can be further constrained. Since the axion field stays in a vacuum manifold 
of U(l) after the spontaneous breaking of the PQ symmetry, the linear topological objects called strings are formed 
in the early universe. Davis (TT| first recognized that axions emitted by these strings give an additional cosmological 
abundance and found that this additional contribution dominates over the coherent oscillations, which gave severer 
upper bound on the axion decay constant. This suggestion is supported by various analytical and numerical studies 
given by different authors jl2Hl4| . However, Harari and Sikivie [15( presented a different scenario which claims that 
the abundance of axions produced by strings do not exceed that of coherent oscillations, and gave a weaker upper 
bound on F a . The subsequent numerical studies provided by [l6j supported this weaker bound. 

This controversy about the contribution from strings arises from different assumptions on the spectrum of radiated 
axions. In refs. flll - flij it is assumed that the typical wavelength of radiated axions is given by the curvature size of 
global strings which is comparable to the size of the horizon, k ~ 27r/£, based on the result of [l7} which claims that 
closed loops or bent stings oscillate many times before they lose most of their energy. Let us call this case as scenario 
A [l8j]. In this scenario, the contribution from axions produced by strings becomes greater than that from coherent 
oscillations by a factor of ln(t/5 s ) w 69, where S s is the width of strings. On the other hand, refs. [H, [l|| suggest 
that the motion and decay of the global strings are more "turbulent". Let us call this as scenario B. In this case, 
strings lose their energy in one oscillation time, quickly decaying into small pieces. Therefore, whole scales between 
the largest scale ~ t and the smallest scale ~ S s give the same contribution to the power spectrum of radiated axions, 
dE/dlnk ~ const, or dE/dk ~ \/k. In this scenario, it turns out that the contribution from strings is comparable 
to that from coherent oscillations. This discrepancy between two scenarios mi ght be resolved by the field-theoretic 
global simulations including the cosmic expansion performed in [l9|, [2(| . In (l9l |2(J| , it was concluded that the power 
spectrum of radiated axions has a sharp peak around the horizon scale fc p hy S ~ 27r/i, supporting scenario A. 

However, it is not sufficient to just consider the string contribution. It was found that these strings are attached 
to surface-like field confi gura tions called domain walls when the axion acquires a mass due to the non-perturbativc 
effect of QCD [2l[ . Lyth |22| pointed out that the annihilation of these domain walls produces additional radiation of 
axions. Subsequently, authors in [23[ and 0] investigated this process, but conclusions of these studies are different 
from each other. Nagasawa and Kawasaki [23| found that axions produced by the collapse of domain walls are mildly 
rclativistic, and this contribution can exceed that from strings. On the other hand, in the study given by Chang, 



Hagmann and Sikivie [24] . the mean energy of axions produced by the decay of domain walls was estimated to be 



larger than that obtained in [23j by a factor of 20. This leads to the conclusion that axions produced by the collapse 
of walls are subdominant compared with that produced by strings. The conclusion of [24| relies on the following 
reasoning. Since domain walls are bounded by strings, the wall energy is converted into the kinetic energy of strings. 
Then, if we assume that scenario B is correct, the spectrum of radiated axions becomes hard (dE/dk ~ l/k). However, 
as we described above, the recent network simulation of global strings supports scenario A. Therefore it is not so clear 
whether the domain wall contribution is significant or not. 
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We point out that this discrepancy on the domain wall contribution is analogous to that on the axionic string 
contribution, and again it can be resolved by preforming the full field-theoretic network simulations. In this paper, we 
aim to determine the contribution of axions produced by the collapse of domain walls bounded by strings and give the 
total relic abundance of cold axions including all production mechanisms (i.e. coherent oscillation, string decay and 
wall decay) . We perform the three-dimensional lattice simulation of the scalar field and follow the whole processes 
relevant to the axion production (from the PQ phase transition to the QCD phase transition). This analysis gives a 
result with least theoretical uncertainties, in the sense that all field configurations arc determined by a first principle. 

We note that the above discussions are only applicable to the case in which PQ symmetry is broken after inflation. If 
PQ symmetry is broken before inflation, there are no contributions from strings and domain walls since the population 
of these defects is diluted in the inflationary era. In this case, isocurvature fluctuations of the axion field gives some 
imprints on anisotropies of cosmic microwave background (CMB) observed today [25|. This observation gives severe 
constraints on axion models and requires significant amounts of fine tunings in the model parameters (called "the 
anthropic axion window") [26j . In this paper, we do not consider this possibility and simply assume that PQ symmetry 
is broken after inflation. 

The organization of this paper is as follows. In section UH we explain the cosmological scenario in which the PQ 
symmetry breaking occurs after inflation and give a qualitative description of topological defects. In section Mil 
we describe the analysis method which wc used in the numerical studies. We present the result of the numerical 
simulations in section IIVI Using the numerical results, we calculate the cosmological abundance of cold axions in 
section [V] Finally, we conclude in section IVT1 



II. COSMOLOGY WITH PECCEI-QUINN FIELD 



In this section, we give an overview of the cosmological aspects of the PQ mechanism, especially concentrating 
on the role of topological defects. Depending on the model parameters, these defects can become either stable or 
unstable. First we introduce theoretical basics and discuss some cosmological consequences. We also give a comment 
on the time scale of the dynamics. 



A. Formation of topological defects and their fates 

We will follow the cosmological evolution of a complex scalar field $ (the Peccei-Quinn field) with the Lagrangian 
density 

C = \\d^\ 2 -V eS [®;T], (1) 

where V^ff[$;T] is the finite-temperature effective potential for the scalar field. At sufficiently high temperature 
(T > F a ), we assume that $ is in the thermal equilibrium, and the effective potential is given by 

V eS [$;T} = ^(\$\ 2 -r, 2 ) 2 + V|$| 2 , (2) 

where we neglect the couplings with other fields for simplicity. Inspection of the form of the effective potential © 
indicates that the PQ phase transition occurs at the temperature 

T = T C = v% (3) 

After that, the scalar field gets vacuum expectation value |(</>)| 2 = rj 2 , and the U(l) symmetry is spontaneously broken. 
Then, due to the spontaneous breaking of the U(l) symmetry, cosmic strings are formed. Because of the causality, 
the population of these strings in the Hubble volume tends to remain in the value of 0(1) (the scaling solution) [271 ]. 
In order to satisfy this scaling property, long strings lose their energy by emitting closed loops of strin gs. These 
loops decay by radiating axion particles with the wavelength comparable to the horizon size [111, [TH, ITa l2dl| . The 
production of axions from decaying string loops continues until the time when the string networks disappear due to 
the mechanism which we describe below. 

When the temperature of the universe becomes comparable to the QCD scale (A ~ lOOMeV), the non-perturbative 
nature of QCD becomes relevant. We can describe this effect by adding the following term in the effective potential ©, 

2 2 

V(8) = ^-(I-cosAdw^), (4) 

JV DW 
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where 9 is the phase direction of the complex scalar field (i.e. $ = |<I>|e ie ), -/Vdw is an integer which is determined 
by considering the color anomaly [2l| . and m a is the mass of the axion. The mass of the axion depends on the 
temperature T, if the temperature is sufficiently high (T > A). Recently, Wantz and Shellard [28| presented the 
temperature dependence of m a which is valid at all temperatures within the interacting instanton liquid model [2£| . 
Fitting the numerical result, they obtained the power- law expression for m a (T) 

A 4 /T\ ~ n 

m a (Tf = c T - - , (5) 



F 2 \A 

where n = 6.68, c T = 1.68 x 1(T 7 and A = 400MeV. This power-law expression should be cut off by hand once it 
exceeds the zero-temperature value m a (T = 0), where 

A 4 

m a {0f=c Q — , (6) 

^ a 

and c = 1.46 x 1CT 3 . 

The existence of the QCD potential (jl]) explicitly breaks the original U(l) PQ symmetry down to its discrete 
subgroup Zjv dw , in which the angular direction possesses the shift symmetry 9 — > 0+27rfc/A r Dw (k = 0, 1, ... , A^dw~ 1). 
This .Zjv D w symmetry is also spontaneously broken because of the vacuum expectation value of the axion field. As a 
consequence, Af DW domain walls attached to strings are formed [2l| . The structure of domain walls depends on the 
number A^dw- If Anw > 1> these domain walls are stable and dominate the energy density of the universe, which 
leads to the discrepancy with the cosmological observations and called the domain wall problem [3(| ■ On the other 
hand, if A^dw = 1, networks of domain walls arc unstable, since the string is attached by only one domain wall. Such 
a piece of the domain wall bounded by string can easily chop the larger one, or shrink itself due to the tension of the 



domain wall [31[ . Hence the networks of domain walls bounded by strings disappear immediately after the formation. 
In this case we can avoid the domain wall problem, and we assume A^dw = 1 hi the rest part of this paper. Note that, 
it is possible to avoid the domain wall problem even in the case with Afow > 1? by introducing the explicit Zn dw 
breaking term in the Lagrangian pH l32j ] . This kind of models also lead to interesting phenomenology [33j , and we 
will present the detailed analysis in another publication [34j . 

In the case with A^dw = 1> the annihilation of string- wall networks occurs around the time of the QCD phase 
transition. At this time, potential energy stored in domain walls and strings is released as radiations of axion 
particles [22|-[24| . As we noted in section |J our interest is to determine whether the population of axions produced 
by this mechanism is comparable or negligible in comparison with other contributions such as axions produced by 
oscillating string loops and the coherent oscillation of the homogeneous field. We will return to this issue in section [V] 
after we present the result of the numerical study in section IPV1 



B. Typical time scale of the dynamics 

Before going to the numerical investigations, let us estimate the typical time scale of the annihilation process. Since 
the tension of domain walls, which causes the decay of string-wall networks, is induced by the existence of the axion 
mass, it is important to note the time at which the axion field begins to "feel" the mass energy. Let us denote this 
time as t\, which is defined by the condition 

m a (Ti) = 3ff(ti), (7) 

where Ti is the temperature at the time t\, and H(ti) is the Hubble parameter at that time. Using the temperature 
dependence of m (T) given by eq. ([5]), we find 

T 1 = 0.981GeV(^) (for T\ > 103MeV), (8) 

or 

Ti = 42.3GeV - " " (^) (for Ti < 103MeV), (9) 

where g„.,i is the radiation degree of freedom at the time t\. Eq. ([S]) is valid only for T\ > 103McV, which corresponds 
to the case in which the condition given by eq. (|7|) is satisfied before m a (T) becomes the zero temperature value m a (0). 
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We must use another expression ((9]) if ?i < 103MeV. However, if we fix the values as g* t i = 70 and A = 400MeV, 
this turnover occurs around the value F a ~ 1.7 x 10 17 GeV. Therefore, we can simply use eq. (J5J as long as we assume 
that F a < 1.7 x 10 17 . The temperature given by eq. ([8]) or cq. ([9]) corresponds to the time 

„/2(4+„)/ Fa \V(4+n)/ A y* 

or 



(l =3.01xlO- s ec(fcL)- (jj^j (_) . (for T, > 103MeV>, (10, 



ii = 1.61 x 10~ 10 scc ( ^ ( ^ . (for T x < 103McV). (11) 

Another relevant time scale is the time when the string-wall networks annihilate themselves. One might guess that 
this occurs when the tension of domain walls dominates over that of strings. Wc denote this time as i 2 , which is 
defined by 

Vwattih) = Mrtr(*2)/*2, (12) 

where cr wa u = 9.23m a (T')F ( 2 is the surface mass density of domain walls HH, /i str = 27rF a 2 ln (^r) is the 

mass energy 

of the strings per unit length, £ is the length parameter of strings defined by eq. ([75)1 . and S s = 1/y/Xrj is the width 
of the core of strings. From eq. (fl"2"j) . we obtain 

, Q „ x -n/2(n+4) / P \ 4 /(™+ 4 ) / A \ ~ 2 

i2 = 1.09xl0-sec(^) (^^) (for T 2 > 103MeV), (13) 

or 

t a = 5.06 x 10- 9 scc f 1Al f°„„ ) ("77^77) (for T 2 < 103MeV), (14) 
and the corresponding temperature 



10 12 GcV/ \400McV 



or 



a ,\-l/("+4) / P \-2/(n+4) / . x 

T 2 = 0.514GeV(^) (^^j (for T 2 > 103MeV), (15) 

where <?* j2 is the radiation degree of freedom at the time £ 2 , and we substituted the typical value In ^^3^) ~ ^9- The 

ratio t 2 /ii ~ 3 indicates that it takes a few Hubble times for the mass term to become effective. 

To summarize, the history of the universe in the case with A^dw = 1 is described as follows. First, the PQ 
phase transition occurs at T ~ T c , and global strings are formed. Next, QCD effects become relevant at T = T\, 
and domain walls are attached to strings. Finally, these domain walls bounded by strings disappear around the 
temperature T = T 2 . Wc will sec that these processes actually occur in the field-theoretic lattice simulations in 
section HVl 



III. ANALYSIS METHOD 



In this section, wc describe the method to calculate the spectrum of axions produced by the decay of string-wall 
networks. Our aim is to extract the pure component of the axion field produced by collapse of the networks, from 
simulated data of the scalar field <£>. In general, the data of 4> contain other components, which can be enumerated 
as follows. 

1. Initial fluctuations. In numerical simulations, we give the initial conditions as Gaussian random fluctuations 
[see eqs. (f3"T)) . (f3"5)) and (|39[) ]. These fluctuations are diluted away by the cosmic expansion, but might not 
be completely negligible even at the final time of the simulation, since the dynamical range of the numerical 
simulation is short. Therefore, they can contaminate the final form of the spectrum of radiated axions. 
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2. Radiations from strings. As we mentioned in section Til Al oscillating loops of strings radiate axions during the 
time between the string formation (T = T c ) and the domain wall formation (T = T±). This contribution must 
be distinguished from the wall-decay contribution which is produced after the time t\ . 

3. Core of defects. In the core of strings, the energy density of the scalar field is higher than that of free axions. 
This can be regarded as another contamination on the spectrum of radiated axions (20| . 

Figure [T] shows the pipeline of removing these contaminations. To remove the contaminations from the core of 
strings, we mask the region near the position of the core of strings and estimate the power spectrum which contains 
only the contribution from free radiations. We calculate the power spectrum in two time slices, the time at which the 
mass of the axion becomes relevant (t = t\) and the time at which the decay of string- wall networks completes (t = to). 
Then, we subtract the spectrum evaluated at t\ from that evaluated at td in order to remove the contributions which 
come from initial fluctuations and radiations from strings. We will give a more detailed description of these procedures 
in the following subsections. First, we establish the formulations for field-theoretic simulations and the notations to 
present the result of the numerical study in section 1111 A[ Then, we describe how to calculate the power spectrum of 
radiated axions in sections IIII Bl and IIII CI Finally, in section IIIIDI wc comment on the subtraction of other radiation 
components. 



t 



— ti initial time of the simulation 




string identification 



masked map 



shift P(fc,ti) 



-ti 



evaluate P 




into 





string identification 
{if it exists) 



and subtract it 






FIG. 1: Schematics of the procedure to estimate the power spectrum of radiated axions. 



A. Formulations 

We assume the flat Friedmann-Robertson- Walker background in which the line element is given by 

ds 2 = dt 2 - R 2 (t)S ij dx l dx j , 



(17) 
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where R(t) is the scale factor of the universe. The equation of motion for <3? is obtained by varying the Lagrangian 
given by eq. (JTJ) with the effective potential given by 

V eS [*;T) = ^(|d>| 2 - ?? 2 ) 2 + ^T 2 |$| 2 +m a (T)W (l - y cos0 ) ' ( 18 ) 

Note that the last term of eq. (|18[) is different from the QCD potential (j4}. We find that the simulation becomes 
unstable when using the potential given by eq. ((H) since this potential is not well-defined at |<f>| = 0. The modified 
potential given by eq. (|18[) avoids this singularity since there is a factor |$| in front of the cosine term. The difference 
between eq. @ and eq. (|18[) is not important in the bulk region on which |$| = 77, and we observe that the quantitative 
behavior of topological defects such as time evolution of the length of strings is unchanged, except the existence of 
the numerical instability. 

We decompose the complex scalar field into its real and imaginary part, such that 

$ = 0i + ifc, (19) 

where cf>i and $2 are real variables. The equations of motion for two real scalar fields cf>i and <p2 are given by, 

& + - |^ = - A<M|$| 2 - V 2 ) - ^T 2 ^ + m 2 (T)ry, (20) 

fa + liHfo - = - A0 2 (|$| 2 - rf) - \t 2 $2- (21) 

If we use the conformal time dr = dt/R rather than the cosmic time, these equations can be rewritten in the following 
form, 

4>'{ - V 2 0x = - A<M|#| 2 - RW) - ^R 2 T 2 ^ + R 3 m 2 a (T)r), (22) 

4% - v 2 2 = - \M\®\ 2 - R 2 v 2 ) - ^R 2 T 2 fo, (23) 

where the prime represents a derivative with respect to r, and we introduced rescaled variable 

$ = (24) 

Note that, in deriving cqs. (|22[) and (|2"3")l . we assume a radiation dominated background, in which R" = 0. 
In the radiation dominated universe, the time and temperature are related by the Friedmann equation 

1 = H2 = SnG^ 4 (5) 
4< 2 3 30 J ' v ' 

where G is the Newton's gravitational constant and is the relativistic degree of freedom. For a convenience in the 
numerical study, we introduce a dimensionless quantity 



45 1 
167r 3 G<7* 77 

Using this parameter, eq. (|25[) can be written as 



C^\/77Z3^-- ( 26 ) 



t=% (27) 

In the following, we normalize all the dimensionful quantities in the unit of r c , which is the conformal time at which 
PQ phase transition occurs [cf. eq. ((3J)], 

$ -> $r c , T -> Tt c , x -> x/t c , etc. (28) 

Also, we introduce the normalized initial Hubble parameter as an input parameter of the numerical simulation 

H(t = U) -> r c H(t = ti) = a, (29) 

and we set the scaling parameter at the initial time into unity 

R(U) = 1. (30) 
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Note that the following relations 



2 C 

Rij c ) = r c /Ti = a and r c = - — . (31) 

6r)a 



By using eq. (|31j) . we can enumerate the various relations in the unit ([28 



l/a, R(t) = or, T 4 = -^L and r? = (32) 

V3 3a 



Then, eqs. of motion for scalar fields (|2"2"]l and (|23p reduce to 



7 2I ,7 h *ia - s \ 4A ^2T , 8r 3 C 3 ^m a (T)^ 2 



$|2 _ 


4r 2 C 2 
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4r 2 C 2 
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V 2 2 = - A0 2 (j$| 2 - j - ^C 2 02. (34) 

Here, the ratio between the axion mass ([5]) and the axion decay constant can be written as 

m a {Tf/ V 2 = c T K n+i (L\ (35) 

where k is the ratio between the QCD scale and the PQ scale, 

k = A/F a = A/??. (36) 

Assuming that $ is in the thermal equilibrium at the high temperature, we give the initial conditions such that the 
two real scalar fields satisfy the renormalized correlation function 

<0 a (x)0 h (y)) = S ab f ^L!*e*<*-y\ (37) 

/d 3 k 
JL^^rik^-^, (38) 

fa.(x)&(y)> = 0, (39) 

where 



1lk = e u k / Ti _ 1 ' Uk = V ^ + m 'ff ' (40) 

and m 2 ff = d 2 V fi/d$*d$>\$ = o is the effective mass of the scalar fields at the initial time. Note that, in eqs. (|37|) and 
(|38p . we subtracted the vacuum fluctuations which contribute as a divergent term when we perform the integral of k. 
In the momentum space, these correlation functions can be written as 

(0 a (k)0 b (k')) = — (2^) 3 ^(k + k'), (41) 

(^(k)^(k')) =<WW2^ 3 5( 3 )(k + k'), (42) 

where </> a (k) is the Fourier transform of 0(x). Since </> a (k) and </> a (k) are uncorrelated in the momentum space, we 
generate </> a (k) and </> a (k) in the momentum space randomly following the Gaussian distribution with 

(l^(k)| 2 ) = ^, (|<Mk)| 2 > = (43) 

and (0 a (k)) = (</> a (k)) = 0, (44) 

for each a = 1 and 2. Then we transform them into the configuration space and obtain the initial field configurations 
</>(x) and </>(x). Here we used (2tt) 3 J( 3 )(0) ~ V, where V = L 3 is the comoving volume of the simulation box and L 
is the size of the simulation box (in the unit of r c ). 

We solve the classical equations of motion given by eqs. ([53"]) and (jM)) in the three-dimensional lattice with 512 3 
points. We impose the periodic boundary condition in the simulation box. The lattice code which we use in this work 
is developed by combining the numerical codes used in [2(| and [33| . We use the forth order symplectic integration 
scheme [36[ to solve the time evolution of the fields. The spatial derivative of the fields is evaluated by using the forth 
order finite difference method. 



B. Energy spectrum of axions 

We calculate the power spectrum of axion radiations P(k, t) defined by 

i<d(i,k)*d(i,k')) = i^(3)(k-k')P(M), (45) 

where (...) represents an ensemble average and d(t, k) is the Fourier component of the time derivative of the axion 
field 

6(t,k) = J d 3 xe* kx d(i,x). (46) 



The value of a(t, x) can be obtained from the simulated data of $ and <j> 



a(t, x) = Im 



The averaged kinetic energy of axions can be written as 



(47) 



On the other hand, the total energy density of axions is given by 

Pa,tot(t) = Pa, kin (t- J ~r Pa,grad(t-J H~~ Pa, mass (t), (49) 

where p Q ,grad(£) is the averaged gradient energy of axions and p a ,mass(t) is the averaged mass energy of axions 



Pa,gr a d(t) = <^|Va(i,x)| 2 ^ , p a , ma ,ss{t) = (^m 2 a a(t 7 x) 2 ^ 



(50) 



One can easily show that, if a(t, x) is a free field, 

PaMn(t) = /0 o ,grad(*) + Pa,mass(i). (51) 

Therefore, P{k, t) can be regarded as the energy spectrum of axions 

/dk 
^P(k,t). (52) 

C. Pseudo power spectrum estimator 

If strings exist, the data of a(t, x) obtained by numerical simulations contain field values around moving strings, 

d(i,x) = df ree (i,x) + (contamination from strings), (53) 

where df ree (£,x) is the contribution from free axion radiations. This moving string contribution can contaminate 
the spectrum of the axion radiations. In order to subtract the contamination from strings we use the pseudo power 
spectrum estimator (PPSE) introduced in [20j | . 

We mask the contribution from the axion field near strings by introducing a window function 

u/{ \ — / (near strings) 

^ ' \ 1 (elsewhere) ' ^ ' 

Then, we obtain the masked axion field 

d(x) = W(x)d(x) = VF(x)d frcc (x), (55) 
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or, in the Fourier space, 

(2tt 

We can compute the power spectrum by using the masked field in a simulation box 



<H k ) = / 7^33 W ( k - k ')a(k')- ( 56 ) 



where V is the comoving volume of the simulation box and flk is a unit vector representing the direction of k. However, 
this masked spectrum is not equivalent to the spectrum of radiated axions, (P{k)) 7^ Ptoc where Pf ree (A;) is defined 
by 

i(a froo (t,k)*a froo (t,k')) = ^^(5 (3) (k-k')Pf r cc(fc,0- (58) 

We can resolve this discrepancy by introducing a window weight matrix, 

1 f dfl k dVLy , 2 

™\"»»)=y2 J ~£ 

and defining the PPSE of P bee (k), 



with M k') satisfying 



M(Kk')^- -±-^\W(k-V)\>, (59) 



k 2 P rik' 

PppsE(fc) = y ^r'ffc, k')P(k% (60) 



t' 2 ^' 9-7T 2 

'l-^-M- 1 (k,k')M(k',k") = ^r5{k-k"). (61) 



Then, we see that (Pp P s E (fc)) = Pfrcc(fc) |2fl|- 

In the numerical simulations, first we calculate the masked power spectrum P(k) and the matrix M[k, k') by using 
the data of $(x) and $(x), then, we compute the power spectrum of free axions by using eq. (1601) . 

For the identification of strings, we use the method developed by 20]. At each lattice segment surrounded by four 
neighboring grids in the simulation box, we identify the existence of the string by the condition A8 > 71*, where A9 
is the minimal phase interval occupied by four grid points in the field space. We compute the position of string as a 
cross-point of two lines with <fii = and cf>2 = in the quadrate. See [20J for details. 

D. Subtraction of pre-existing radiations 

In order to extract the spectrum of axions radiated from the decay of string-wall networks, we must subtract the 
contributions of pre-existing radiations, which have been created before the decay of domain walls, from the whole 
spectrum calculated by using the field data obtained after the decay of networks. The radiations created before 
the decay of domain walls are diluted due to the cosmic expansion, and we evaluate this redshift factor before we 
perform the subtraction. Note that the axion mass becomes non-negligible around the time of the decay of domain 
walls. Hence, we cannot subtract the spectrum simply assuming that the spectrum is diluted as R~ 4 , which is only 
applicable to massless particles. 

From eq. (|52|). the total energy density of axions can be written as 

d 3 k 

where 



A»(*)= / T^rgMMWM), (62) 



w Q (M) = ^m 2 a + k 2 /R{tf (63) 
is the energy of axions with momentum k/ R(t), and we define 



(64) 
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We can regard ri a (k,t)d 3 'k/ (2tt) 3 as the number density of axions which have comoving momentum within the range 
from k to k + dk. Therefore, we expect that n a (k, t) scales as R(t)~ 3 , if there are no changes in the number of axions. 

Let us denote the time at which domain walls sufficiently decay as td (this definition of td may contain an ambiguity, 
which we discuss in the next section). By using the fact that n a {k, t) oc R(t)~ 3 if there is no absorption or production 
of axions, and eq. (|64|) . we find the form of the spectrum of pre-existing radiations produced before t\, where t\ is 
defined by eq. {7J, at the time td, 

P rad(Md ) = mh)^^(^) 3 , (65) 

where P(k,t\) is the spectrum evaluated at t\. Subtracting the contribution P ra d{k,td) from the whole spectrum 
P(k,td) evaluated at td, we obtain the spectrum of radiations produced after the decay of domain walls 

Pdcc(Mrf) = P(k,U) - P xad (k,t d ). (66) 



IV. RESULTS OF NUMERICAL SIMULATIONS 



In numerical simulations, we can vary four parameters A, n, C an d ol. We set a = 2.0 which corresponds to the fact 
that Tj = 0.5t c . Also, we choose the value of ( as 3.0, which corresponds to the conditions r\ = 1.23 x 10 17 GeV and 
g* = 100. It seems that this value of r\ may be too high and affect the small scale dynamics of the system. We will 
discuss this point in the end of this section. Note that, from eq. (|3"2"f we see that 77 = 1 in the unit of t" 1 . Other 
parameters that we used are summarized in table H] The dynamical range of the simulation is estimated as Tf/ri = 24. 



TABLE I: Parameters used in numerical simulations. 



grid size (TV) 


512 


box size (L) 


20 


total number of steps 


1150 


time interval (dr) 


0.01 


A 


1.0 


K [eq. 


varyin] 


C [eq. (HJlj 


3.0 


a [eq. JUJ] 


2.0 


c T [eq. ©j 


6.26 


c [eq. ©] 


1.0 


initial time (n) 


0.5 


final time (77) 


12.0 



We must keep the following conditions, in order to simulate the dynamics of the topological defects correctly. 

• The width of global strings S s = l/rjy/X must be greater than the physical lattice spacing Sx p hys = R(t)L/N, 
where N = 512 is the number of grids, in order to maintain the resolution of the width of strings. 

• The Hubble radius H -1 must be smaller than the box size R(t)L, to avoid the unphysical effect caused by the 
finitcness nature of the simulation box. 

The physical scale of the Hubble radius and the width of strings divided by the physical lattice spacing are respectively 

H- 1 N 5s 3iV 

7 = -r T i and 7 = t=- (0 7 ) 

ozphys L ox phys 2LC r vA 

For the parameters given in table HI we get H~ l /5x v \ lys ~ 307 and <5 s /fe P hys — 1-07 at the end of the simulation 
t = Tf. Therefore, the conditions described above are satisfied even at the end of the simulation. 

We also treat cr and Co defined in eqs. ([5]) and ^ as free parameters in numerical simulations. In terms of these 
parameters, the time at which the value of m a (T) reaches the zero temperature value m a (0) is written as 

r a = 1.73 x ( -M kT 1 . (68) 
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Here we used the conformal time with the unit of r c = 1. Also, the time t\ given by cq. (JT)) is rewritten as 

Tl = 1.52x(Mj c- 1/{n+A) K-K (69) 

Choosing the values of three parameters cq, ct and K corresponds to the fact that we tune the values of r a , t\ and 
TO a (0) in the numerical simulations. Unfortunately, due to the limitation of the dynamical range, we cannot choose 
all parameters to be realistic values. One possible choice is to fix the ratio between t\ and r a so that 

/ \l/n / „„ N 2/(n+4) /, n \2/(«+4) 

n/r. = 0.88 x r^j c?'M(?f) =0 - 97x (t) ' (70) 

where the last equality follows from the realistic values ct = 1-68 x 10 -7 and Co = 1.46 x 10~ 3 . In order to satisfy 
the condition (|70p . we must choose ct and cq so that 

c T c - (n+4)/4 ~ 6.26. (71) 

In the numerical simulations, we use the values Co = 1.0 and ct = 6.26 which satisfy the above condition. However, 
we note that in this case the expression for t-i is given by eq. (|14l) rather than eq. (|13|) . 

r 2 = 1-65 x(^) (£) V2 c^'\-\ (72) 

since the value of T2 turns out to be greater than r . In other words, in the numerical simulations the condition given 
by eq. (|12[) is satisfied after the axion mass m a (T) reaches the zero-temperature value m a (0), while in the realistic 
case with F a < 1.7 x 10 17 it occurs before m a {T) reaches m o (0). This discrepancy is caused by the smallness of (. 
Namely, the realistic value of £ ~ Mpjr\ is much greater than 3.0 which we used in the numerical simulations, where 
Mp is the Planck mass. The expression (jT2"j) also depends on /3 = ln(t/\/f<5 s ) which comes from /i str in eq. (fTS")) . Here 
we use the value (3 ~ 4 obtained by using parameters which are used in numerical simulations. We summarize the 
typical time scales given by eqs. ([68]), P5) and $2$ in table El 



TABLE II: Typical time scales for various values of k. 





Tl 




T2 


0.4 


3.20 


3.29 


4.12 


0.35 


3.65 


3.76 


4.71 


0.3 


4.27 


4.38 


5.50 


0.25 


5.12 


5.26 


6.60 


0.2 


6.40 


6.57 


8.25 



Now, let us show the results of the simulations. Figure [2] shows the visualization of one realization of the simulation. 
We see that at the first stage of the simulation, strings evolve with keeping the scaling property However, at late 
time they shrink because of the tension of domain walls. We also show the spatial distribution of the phase of the 
scalar field 9 in figure[3] Note that the width of domain walls ~ m^ 1 is much greater than that of strings ~ (v^A??) -1 , 
as shown in figure [3] 

We performed 20 realizations for each choice of the parameter n. For each realization, we calculated the length 
parameter of strings 

i = P^l t ^ (73) 

Mstring 

and the area parameter of domain walls 

A = — t. (74) 

Cwall 

Figure @] shows the time evolutions of £ and A for various values of k. Comparing the plot of A with tabic UH we see 
that the value of A deviates from the scaling behavior (A ^constant) around r = n and begins to fall off around 
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t = T2- Note that £ starts to fall later than A docs. This can be interpreted as follows. Since domain walls are 
two-dimensional objects, they curve in various directions. This curvature gets stretched when the tension of walls 
becomes effective. The stretching process of walls reduces the value of A, but might not affect the length of strings 
(i.e. the value of £). Later, stretched walls pull the strings attached on their boundaries, which causes the reduction 
of£. 
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FIG. 2: Visualization of one realization of the simulation. In this figure, we take the box size as L = 15 and N = 256, which is 
smaller than that shown in table [I] Other parameters are fixed so that A = 1.0, £ = 3.0, a = 2.0 and k = 0.4. The white lines 
correspond to the position of strings, while the blue surfaces correspond to the position of the center of domain walls. 

We also calculated the spectrum of axions radiated from strings and domain walls, using the method described in 
previous section. Figure [5] shows the spectra of free axions evaluated at t\ and td- The basic behavior of the spectrum 
evaluated at t\ is similar to that obtained in ref. (20j . This spectrum is dominated by the contribution of axions 
produced by strings. However, the population of axions with high momenta increases after the decay of domain walls 
(t = td)- The final form of the spectrum, obtained by subtracting the components of radiations produced before ti, is 
shown in figure O The spectrum has a peak at the low momentum. This disagrees with the result of Chang, Hagmann 



14 





FIG. 3: The distribution of the phase of the scalar field 9 on the two-dimensional slice of the simulation box. In this figure, 
we used the same data that are used to visualize the result with r = 5.0 in figure [21 The value of 9 varies from —it (blue) to n 
(red). Domain walls are located around the region on which 9 passes through the value ±ir, while the green region corresponds 
to the true vacuum (9 = 0). The length scale of the change of 9 is roughly estimated as ~ m" 1 , which gives the thickness of 
domain walls. 
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FIG. 4: Time evolution of the length parameter £ (left panel) and the area parameter „4 (right panel) for various values of n. 
Although walls do not exist before the time n, we can show the value of A evaluated at the time r < Ti. This is because the 
value of A is calculated from the number of grid points on which the phase of the scalar filed passes the value 9 = tt. In this 
sense, A represents the area of domain walls only after the time ti. 



and Sikivie [24[, which claims that the radiated axions have a spectrum proportional to 1/k. Note that, however, 



there is a high frequency tail in the spectrum, which has a cutoff at the momentum corresponding to (twice the size 
of) the width of strings k ~ (2ir/25 s )R(t ( i) — 64.4 (for k = 0.3). This feature might be interpreted in terms of the 
reasoning of [HI, HUE! (scenario B). Namely, there are various size of defects around the time t = td, and small scale 
defects can radiate axions with harder momentums. As we see in figure |H1 the contribution from these hard axions is 
subdominant, and most axions have a momentum comparable to the mass of the axion k/R(td) ~ m a . 



15 



Using the result of Pd ec (k-,td)-, we compute the mean comoving momentum of radiated axions 

J§P dec (k,t d ) 



k(t B 



We define the ratio of the physical momentum kj R{td) to the axion mass m a (td), 

= ~k{td)/R{t d ) 
m a (t d ) 

From the result of the numerical simulations with k = 0.3, we obtain 

k(t d ) = 5.76 ±0.15 and e w = 3.12 ± 0.08. 
This result corresponds to the mean energy of axions [cf. eq. ([63p], 



C0a(t d )/m a {t d ) = y/l + el = 3.28 ± 0.08. 



(75) 
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FIG. 5: The power spectrum of free axions calculated by using PPSE formula (|60[l in the simulations with k = 0.3. We plot the 
spectra evaluated at two different times t\ and td- Note that the result of P(k,ti) shown here does not contain the numerical 
factor defined in eq. ([65)l 

Note that there are some ambiguities in this analysis. For instance, we choose the time td, at which the decay of 
networks completes, by hand. If we choose td as sufficiently late time (for example, the final time of the simulations), 
we would underestimate the mean momentum of radiated axions defined by eqs. (|75[) and (|76|). since the physical 
momentum gets rcdshiftcd proportionally to 1/R(r d ). Figure [7] shows the results of the physical momentum k/R(rd) 
for various choices of Td in the simulations with k = 0.3. The value of k/R(rd) begins to shift as oc l/R at the result 
with Td = 10.25. This value of T d corresponds to the time at which the area of domain walls becomes less than 0(1)% 
of the Hubble scale (A <0.01). Therefore, we choose Td as the time at which the value of A falls below 0.01. We use 
this criterion to calculate the spectra with different values of n shown in figure [5] 

Another subtlety is whether the results of numerical simulations are sensitive to the choice of k = A/F a . There is a 
tremendous hierarchy between the QCD scale and the PQ scale, A/F a ~ 100McV/10 10 GcV = 10 -11 , but we cannot 
perform the simulations with such a small value of k because of the limitation of the dynamical range. Nonetheless, 
we believe that the ratio between the mean momentum of the radiated axions and the typical momentum scale (such 
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FIG. 6: The spectrum of axions produced by the decay of networks, defined by eq. (|66l) . Note that the results with different 
value of k are evaluated at different times (jd). The form of the spectra is different from the relation Pdec(k) oc 1/k which is 
indicated by the dotted line. 
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FIG. 7: The value of the mean physical momentum of radiated axions k/R{jd) and its dependence on the choice of r d . The 
dashed line represents the relation k/R(rd) oc 1/R(r d ) oc l/r d . The five points correspond to the results with the value of the 
area parameters (a) A{r d ) <0.1, (b) A{r d ) <0.05, (c) A(r d ) <0.01, (d) A{r d ) <0.005 and (e) A{r d ) <0.001. 
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as m a ) is not so sensitive to the value of k, since the power spectrum has a sharp peak at the typical momentum scale 
as we see in figure [51 However, one might regard the peak momentum as the inverse of the horizon scale ~ 2n/td, 
instead of the mass of the axion. The wrong interpretation of the peak momentum affects the estimation of the present 
abundance of axions, which will be given in the next section. We present the result of the ratios [k/R(td)]/m a (td) 
and \k/ R{td)]/ /td) for k = 0.3, 0.35 and 0.4 in table IIII1 We observe that the K dependence of [k/R(td)]/m a {td) is 
weaker than that of [k / R(td)] / (2tt / td) ■ This fact implies that it is reasonable to assume k/R(td) oc m a {td) rather than 
k/R(td) oc 2-rr/td- For now, it is difficult to discuss the significance of the result in which the value of [k / R(td)] /m a {td) 
varies with «, because of the lack of samples. It can be said that the result with larger value of k is unreliable 
since the networks of strings begin to collapse soon after their formation. Indeed, in figure |6j we observe that the 
high- momentum cutoff, which is the consequence of the "smallness" of k or m a compared with F a , is less apparent 
in the result with k = 0.4. This ambiguity should be resolved in future numerical studies with improved dynamical 
ranges. Regarding this ambiguity, we use the result with the smallest value of k given by eq. ([77)) or (|78[) . 

TABLE III: Ratio between the mean momentum of radiated axions k/R(t d ) and the mass of the axion m a (td), or the mean 
momentum of radiated axions k/R(td) and the inverse of the horizon scale 2ir/t d , for different values of 

k [k/R(t d )]/m a {t d ) {k/R(td)}/(2ir/t d ) 

0.3 3.12±0.08 4.70±0.12 

0.35 2.64±0.09 3.71±0.13 

0.4 2.51±0.10 3.13±0.13 



Furthermore, there are ambiguities in the values of scaling parameters defined in eqs. (|73j) and (1741). Our result 
£ ~ 0.5, shown in figure^] is somewhat lower than the value £ ~ 0.8-1.3 obtained in previous studies [19l. I20I l37l [38j . 
This might be caused by the different choice of the parameter £ used as an input of the numerical simulations. Our 
choice ( = 3.0 is smaller than the values C = 8-10 used in past numerical simulations [H, H3, IH[. The parameter £ 
controls the magnitude of the symmetry breaking scale 77 [see eq. (|32[) ]. which determines the width of global stings 
S s oc I/77. Therefore, different choice of £ affects the emission rate of Nambu-Goldstone bosons from strings (27L l39l| 
Tng = r/[27rL s ln(L s /(5 s )], where T is a numerical factor of 0(10-100) and L s ~ t is the characteristic length scale 
of strings. The simulation with small value of C corresponds to the simulation with thick strings, in which the global 
string networks lose their energy efficiently due to the emission of Nambu-Goldstone bosons, since the logarithmic 
correction to the emission rate T^g oc 1/ h\(t/S s ) becomes large. This large emission rate of Nambu-Goldstone bosons 
reduces the energy density of global string networks and suppresses the value of £. However, it was argued that in the 
realistic case with ln(t/d s ) « 70 the radiative effect becomes subdominant, and the value of £ is purely determined by 
the formation rate of loops [H, Hfj. Regarding this effect, the authors of ref. [38| estimated the final value of scaling 
parameter as £ = 1.6±0.3. Indeed, the results with smaller values of k in figure [4] indicate that the value of £ increases 
due to the change of emission rate Tng oc 1/ ln(t/5 s ) with time. We anticipate that the value of the length parameter 
gradually reaches the final value £ ~ 1, which cannot be observed in the simulations with the limited dynamical range. 

The reason why we choose smaller value of £ is to improve the dynamical range of simulations by keeping the width 
of strings greater than the lattice spacing [see eq. ([57)1 ]. This choice enables us to try to perform simulations with 
varying the values of k but invalidates the estimation of £ due to the large emission rate of Nambu-Goldstone bosons. 
However, we believe that this choice does not affect our main result about the radiated axions produced by domain 
walls, since the choice of the parameter £ only controls the small scale properties such as the width of strings, while 
the population of axions are dominated by low momentum modes which produced by the large scale physics with the 
wavelength comparable to the inverse of the axion mass. 

The precise determination of the values of scaling parameters including the effect of backreaction of Numbu- 
Goldstone boson emission is beyond the scope of this paper. To be conservative, we use the rough estimate £ ~ 1.0±0.5 
with 50% uncertainty when we calculate the abundance of cold axions in the next section. We also assume that the 
area parameter A possesses similar uncertainty, and use the value A ~ 0.50 ± 0.25 around the time of the formation 
of domain walls. 



V. RELIC ABUNDANCE OF COLD DARK MATTER AXIONS 

In this section, we calculate the abundance of cold dark matter axions. These axions are produced in three processes: 
(A) vacuum misalignment, (B) string decay and (C) domain wall decay In the following, we treat these three cases 
respectively and estimate the energy density of axions at the present time. 
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A. Zero modes 



The averaged homogeneous value of the axion field (zero mode) begins to oscillate around the minimum of the 
potential at the time t\. Let us denote the initial amplitude of the 9 angle as 6\. The energy density of the zero 
modes is given by 

Pafi{tx) = (79) 

Noting that the number of zero mode axions is fixed after the time t\ , we find the energy density at the present time 
to 

Pa,o{to) = Pa,o{ti) — FFT ^77T • ( 80 > 
m a (Ti) \R(t Q ) J 



From the entropy conservation, it follows that 

/i?(ti)\ 3 s 



where sq is the entropy density at the present time, which satisfies 

sgh 2 _ 4 g*s.o M R h 2 
Pc.o 3 g*,o T 



(81) 



(82) 



Here, p c .o is the critical density today, g*s,o and <?*,o are the effective degrees of freedom for entropy density and 
energy density of radiations at the present time [4l[, T is the temperature today, flnh 2 = pFt(to)h 2 / p c ,o is the density 
parameter of radiations, and h is the reduced Hubble parameter: H — 100/ikm-scc _1 Mpc _1 . Using eqs. (|79")l - ([52")) 
and the expression for T\ given by cq. ([5]), we find that the density parameter of the zero mode axion n a Q h 2 = 
Pa,o(to)h 2 1 p cfi becomes 

n a0 h 2 = 0.095 x 9 2 f^i) — . 83 

If PQ symmetry is broken after inflation, the value of 9\ spatially varies, and we can replace 9\ by the root-mean- 
square value 

{ ° 2i) = hf 02id01 = n2/3 - 

Furthermore, it was pointed out that the anharmonic effect becomes important for a large value of 0\ (42Tl44l |. 
Considering this effect, we take the replacement 6\ — > Canh^g-j where c an h is the anharmonic correction. Turner [42[ 
calculated the anharmonic effect numerically and obtained the correction factor c an h = 1.9-2.4. Later, Lyth [43| gave 
the extensive calculation and reported the agreement with Turner's result within a factor of 2. The more precise 
calculation given by [44[ leads c an h — 1-85. Here we take the value Q\ — > 1.85 x and obtain 

(n+2)/2(n+4) / F a \ ( n+6 )/(™+ 4 ) / A \ 



fl a0 h 2 = 0.58 x ( ^i) — . 84 

V 70 / ^10 12 GeVy V 4 00MeVy v ; 



B. Axions radiated from scaling global strings 

We can estimate the energy density of axions radiated from strings by using a similar method introduced in appendix 
B of [20| . In [20( 1 . it was assumed that the number of axions is fixed after the time ti. However, since the mass of 
the axion becomes non-negligible after the time t\ (recall that t\ < £2), the analysis of [2(| is applicable only for the 
radiations produced before the time t\. Here we assume that the radiation of axions from scaling global strings is 
terminated at t\. The abundance of axions radiated after t\ will be estimated in the next subsection. 

The present number density of axions radiated before t\ is estimated as [20l | 

fR(h)\ 3 2F 2 f (h/Vt\ ( ^ 

na > M= {-m) —-e ln {—)> (85) 
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where e is the ratio between the mean momentum of axions radiated by strings and the horizon scale. 

The energy density of axions today is given by p a .str(^o) = TO a(0)«. a ,str(^o)- Then, we found the density parameter 
of axions radiated by strings 

£ / 5a<1 N-(n+2)/2(n+4) / Fa \ (n+6)/(n+4) , A v 



In numerical simulations performed by ref. [13], the value of e is estimated as e _1 = 0.23 ± 0.02. By using this value 
for e and the rough estimation for the length parameter £ ~ 1.0 ± 0.5, we obtain 



tt astr h 2 = (4.0 ±2.0) x (^i) ( — ^ ) ( I. (87) 

a ' str v ' V 70 / V 10 12 GeV J \ 400MeV / v ; 



( ra +2)/2(n+4) / ^ a \ ("+6)/(n+4) / A \ 



V400MeV / 



C. Axions radiated from the decay of string- wall systems 



Define the area parameter of domain walls at t = t\ 



and the length parameter of strings at t = t\ 



. _ Pwall(il) 

Ai = 7TT 1 ' 

0"wallitlj 



Pstring(il) ,2 



ti = (89) 



The string- wall networks begin to collapse around the time i = ti. We simply assume that, after the time ti, the 
whole energy stored in these defects is diluted as R(t)~ 3 due to the cosmic expansion 



Pstring— wall (t) 



. CTwall(il) t /^string 

(*i) 

•Ai h ?i — 



for t>h. (90) 



Suppose that the decay is complete at the time td > t%. The number density of axions produced by the decay of 
string- wall networks is 



n 



a, dec (t) 



_ Pstring— wall (td) ( R{td) 



0J n 



\R(t) 



1 



y/1 + elm a (t d ) 



/^string 

(ti) 

A-i h £i 



'1 



(91) 



where Q a = yT + e l] m a{td) is an average of the energy of radiated axions [see eq. ([75])]. 

The above expression does not depend on td except the factor l/m a (td). However, since the change in the mass 
of the axion can be negligible (|m a /m^| ~ H/m a < 1) for t > t\, we can approximate m a (td) ~ m a (ti). Then, the 
present energy density of axions radiated after t\ is given by 



m Q (0) 

Pa,dcc{to) = m a {Q)n a ^ cc (t ) - 



. CTwall(tl) /^string 

(*i) 

Al h ?1 72 



\R(t ) 



(92) 



The density parameter of axions radiated from the decay of defects is given by 

W = 8.46 x 10- x -■«>«»«> /^y«>«-«> 

V / T+1^ V 70 / Vl0 12 GeVy V400MeVy v y 

As we discussed in section HVl we use the conservative estimations £i ~ 1.0 ± 0.5 and A\ ~ 0.50 ± 0.25. Substituting 
these values and the value of e w given by eq. (|77p. we finally obtain 



Comparing eqs. (|87|) and (|94[) . we see that the contribution from domain wall decay is greater than that from string 
decay. This result supports the conclusion of refs. [22[ and (23|. 



(n+2)/2(n+4) / F a \ ("+ 6 )/("+ 4 ) / A \ 
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VI. CONCLUSION 

We have investigated the production of axions from annihilation of domain walls bounded by strings. We followed 
the evolution of topological defects by solving classical field equations on the three-dimensional lattice. In numerical 
simulations, we observed that global strings indeed annihilate when the mass of the axion becomes relevant. We 
calculated the power spectrum of axions produced by the annihilation of string-wall networks by subtracting the 
contribution which contains fluctuations given by initial conditions and radiations produced by oscillating strings. 
The spectrum has a peak at the low frequency, and the mean energy of radiated axions is ui a ~ 3m Q , which is 
consistent with the result of [23| . 

The total abundance of cold dark matter axions is given by the sum of eqs. (|84|) . (|87|) and ([94]), 

^a,tot^ 2 = ^a,0^- 2 + ^a,str^ 2 + ^a,dec^ 2 

= < 16±6 K^f IsO^)- (95) 

where we used n = 6.68 according to [HI]. The large uncertainty arises from the poor determination of the scaling 
parameter £ ~ 1.0 ± 0.5, which might be fixed by developing the model of the evolution of global string networks, 
such as the study given by j3{|. We require that Q a ,toth 2 must not exceed the observed value of the abundance of the 
cold dark matter fi^DM^ 2 = 0.11 0- This gives an upper bound for the axion decay constant 

F a < (1.2-2.3) x 10 10 GeV, (96) 

if we take g*.\ = 70 and A = 400MeV. This bound is severer than the result of the previous study F a < 3x 10 11 GeV 1201 . 
which is obtained by considering only the abundance of axions radiated by strings. We note that the other group [28[ 
already reported another bound F a < 3.2^ x 10 10 GeV as severe as obtained here, although they considered only the 
contribution of axionic strings. We believe that this severity would come from the larger scaling parameter £ « 13 
used in their analysis than our numerical prediction £ ~ 1, which overestimates the relic abundance of axions radiated 
by strings. In this sense, we regard that the axion string constraint is milder than that indicated by eq. (|96|) . 

As we mentioned in section U there is a lower bound F a > 10 9 GcV which comes from astrophysical observations. 
Combining this lower bound with the bound (|96[) . we conclude that axion models are constrained into the narrow 
parameter region F a ~ 10 9 -10 10 GeV, which corresponds to the axion mass m a ~ 10 _3 -10 _2 eV. 

We note that our numerical result contains some ambiguities as we discussed in section IIVI Especially, we choose 
the unrealistic values of parameters such as ft, ct and Co, which determines the magnitude of the axion mass through 
eqs. ([5]) and ©, in order to follow the whole relevant processes within one realization of the numerical simulation. It 
is not obvious whether the results are sensitive to the tuning of these theoretical parameters. This should be tested 
in future high-resolution numerical studies with larger dynamical ranges. 
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